"""This demo program solves Poisson's equation

    div C grad u(x, y) = f(x, y)

on the unit square with source f given by

    f(x, y) = - 2*(1 + pow(x[0],2) + x[0]*x[1] + pow(x[1],2))*exp(x[0]*x[1])

and boundary conditions given by

    u(x, y) = exp(x*y)  for x = 0 or x = 1 or y = 0 or y = 1

The conductivity C is a symmetric 2 x 2 matrix. The conductivity is

    C = ((2.0, 1.0), (1.0, 2.0))

This corresponds to the example 2 of Darlan Thesis (printed copy, p.68)

"""

from dolfin import *

# Read mesh from file and create function space
mesh = UnitSquare(4, 4)
V = FunctionSpace(mesh, "CG", 1)

# Define Dirichlet boundary (x = 0 or x = 1 or y = 0 or y = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or \
           x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Expression('exp(x[0]*x[1])')
bc = DirichletBC(V, u0, boundary)

# Define conductivity matrix
C = as_matrix(((2.0, 1.0), (1.0, 2.0)))

# Define variational problem
v = TestFunction(V)
u = TrialFunction(V)
f = Expression("- 2*(1 + pow(x[0],2) + x[0]*x[1] + pow(x[1],2))*exp(x[0]*x[1])")
a = inner(grad(v), C*grad(u))*dx
L = v*f*dx

# Compute solution
problem = VariationalProblem(a, L, bc)
u = problem.solve()

# Save solution in VTK format
file = File("poisson.pvd")
file << u

# Calculate max error norm
u0_interp = interpolate(u0, V)
u_exact = u.vector().array()
u_numerical = u0_interp.vector().array()
diff = abs(u_numerical - u_exact)
print 'Max Error: %1.2E' % diff.max()

# Plot solution
plot(u, interactive=True)
